library(ggplot2)
library(haven)
library(list)
library(dplyr)
library(stargazer)
library(estimatr)
library(ggpubr)
library(scales)
library(xtable)
library(MASS)
library(ordinal)


#####################################################################################################################
# Replication code for "Corrupted Estimates? Response Bias in Citizen Surveys on Corruption", Political Behavior
# Author: Mattias Agerberg, University of Gothenburg, mattias.agerberg@gu.se
# MAIN MANUSCRIPT
#####################################################################################################################


#####################################################################################################################
# Paper wd
setwd("") 

# Data
data <- read.csv("data_RO.csv") 

#####################################################################################################################
# ADDITIONAL DATA MANAGEMENT
#####################################################################################################################
attach(data)

# Direct bribe question
data$bribe[Bribe_direct=="Nu"] <- 0
data$bribe[Bribe_direct=="Prefer sa nu raspund"] <- 0
data$bribe[Bribe_direct=="Da"] <- 1

# Reverse code Corruption change and Economy - high values: corruption increase, worse economy
data$corr_soc <- data$corr_soc*(-1)+6
data$econ <- data$econ*(-1)+6

# List variables
list <- cbind(List_nonsensitive,List_sensitive)
data$items[is.na(list[,1]==F) | is.na(list[,2]==F)] <- rowSums(list, na.rm = TRUE)
data$items <- as.integer(data$items)

data$list_treatment <- NA
data$list_treatment[is.na(List_nonsensitive)==F] <- 0
data$list_treatment[is.na(List_sensitive)==F] <- 1
data$list_treatment <- as.integer(data$list_treatment)

data_complete <- data[,c("items","list_treatment","gov_support_01","female","city_over200k","university_ed","Age","income_top20")]
data_complete <- data_complete[complete.cases(data_complete),]

#####################################################################################################################
# LIST ANALYSIS
#####################################################################################################################
# Difference-in-means
fit.list <- ictreg(items ~ 1,
                   data = data, treat = "list_treatment", J = 4, method = "lm")

fit.sens <- glm(bribe ~ 1,
                data = data, family = binomial("logit"))

avg.pred.social.desirability <- predict(fit.list,
                                        direct.glm = fit.sens, se.fit = TRUE)

# FIGURE 3
pred.bias <- as.matrix(avg.pred.social.desirability[[1]])
pred.bias <- data.frame(pred.bias)
pred.bias$est <- c("List estimate","Direct estimate","Difference")
pred.bias$est <- factor(pred.bias$est)
pred.bias$est = factor(pred.bias$est,levels(pred.bias$est)[c(2,3,1)])

ggplot(pred.bias,aes(x=est,y=fit,colour=est)) +
  geom_errorbar(aes(ymin=lwr, ymax=upr), width=.1,position = position_dodge(width = 0.2)) +
  geom_point(position = position_dodge(width = 0.2),size=3.5) +
  scale_y_continuous(limits=c(0,0.45)) +
  labs(y="Estimated Proportion",x="") +
  theme_minimal(base_size = 11) + scale_colour_grey(start=0.5, end=0.2) + 
  theme(legend.position="none", axis.title.y = element_text(vjust = 2.4),axis.text.x = element_text(size=10.5))
ggsave("listbias.pdf")

# TABLE 2
list.table <- matrix(nrow = 3,ncol = 3)
colnames(list.table) <- c("Proportion","95% CI low","95% CI high")
rownames(list.table) <- c("List estimate","Direct estimate","Difference")
for (i in 1:3) {
  list.table[1:3,i] <- as.numeric(pred.bias[1:3,i])
}

# Tablenotes were added using "threeparttable" in LaTeX
stargazer(list.table, title = "Sensitive item (being asked to pay a bribe), estimates", rownames = TRUE, align = TRUE,
          notes = "`List estimate' indicates the main estimate from the list experiment, obtained via the difference-in-means estimator. `Direct estimate' is the result obtained for the direct question using a logistic regression model, and `Difference' is the difference between the two estimates.")

# FIGURE 4
# List analysis, subgroups fit
fit.list <- ictreg(items ~ gov_support_01 + female + city_over200k + university_ed + income_top20 + Age,
                   data = data_complete, treat = "list_treatment", J = 4, method = "nls")

fit.sens <- glm(bribe ~ gov_support_01 + female + city_over200k + university_ed + income_top20 + Age,
                data = data, family = binomial("logit"))

# List analysis, Gender
data_complete_female <- data_complete_male <- data_complete
data_complete_male[, "female"] <- 0
data_complete_female[, "female"] <- 1

avg.pred.male <- predict(fit.list,
                              newdata = data_complete_male, avg = TRUE, interval = "confidence")

avg.pred.female <- predict(fit.list,
                                 newdata = data_complete_female, avg = TRUE, interval = "confidence")

newdata2 <- with(data, data.frame(female = 1, gov_support_01 = mean(gov_support_01), city_over200k = mean(city_over200k), 
                                  university_ed = mean(university_ed), income_top20 = mean(income_top20), Age = mean(Age)))
pred.sens.female <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

newdata2 <- with(data, data.frame(female = 0, gov_support_01 = mean(gov_support_01), city_over200k = mean(city_over200k), 
                                  university_ed = mean(university_ed), income_top20 = mean(income_top20), Age = mean(Age)))
pred.sens.male <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

# Plot by gender
pred.mf <- matrix(nrow = 4,ncol = 5)
pred.mf[c(1,3),4] <- "Female"
pred.mf[c(2,4),4] <- "Male"
pred.mf[1:2,5] <- "List"
pred.mf[3:4,5] <- "Direct"
for (i in 1:3) {
  pred.mf[1,i] <- as.numeric(avg.pred.female$fit[[i]])
}
for (i in 1:3) {
  pred.mf[2,i] <- as.numeric(avg.pred.male$fit[[i]])
}
pred.mf[3,1] <- pred.sens.female[[1]]
pred.mf[3,2] <- pred.sens.female[[1]]-1.96*pred.sens.female[[2]]
pred.mf[3,3] <- pred.sens.female[[1]]+1.96*pred.sens.female[[2]]
pred.mf[4,1] <- pred.sens.male[[1]]
pred.mf[4,2] <- pred.sens.male[[1]]-1.96*pred.sens.male[[2]]
pred.mf[4,3] <- pred.sens.male[[1]]+1.96*pred.sens.male[[2]]

pred.mf <- data.frame(pred.mf)
pred.mf$X1 <- as.numeric(levels(pred.mf$X1))[pred.mf$X1]
pred.mf$X2 <- as.numeric(levels(pred.mf$X2))[pred.mf$X2]
pred.mf$X3 <- as.numeric(levels(pred.mf$X3))[pred.mf$X3]

gender.list.gr <- ggplot(pred.mf,aes(x=X4,y=X1,colour=X5)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.25) +
  geom_point(position = position_dodge(width = 0.2),size=1.3) +
  scale_y_continuous(limits=c(0, 0.83)) +
  labs(y="",x="",colour="") +
  theme_minimal(base_size = 12) + scale_colour_grey(start=0.5, end=0.2) +
  theme(axis.text.x = element_text(size=6.5), axis.title.y = element_text(size=5.5, vjust = 1), 
        axis.text.y = element_text(size=5.3), legend.text = element_text(size=6.5),
        panel.grid.minor = element_line(size = 0.15), panel.grid.major = element_line(size = 0.3))

# List analysis, Gov supporters
data_complete_nongov <- data_complete_gov <- data_complete
data_complete_gov[, "gov_support_01"] <- 1
data_complete_nongov[, "gov_support_01"] <- 0

avg.pred.gov <- predict(fit.list,
                         newdata = data_complete_gov, avg = TRUE, interval = "confidence")

avg.pred.nongov <- predict(fit.list,
                           newdata = data_complete_nongov, avg = TRUE, interval = "confidence")

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = 1, city_over200k = mean(city_over200k), 
                                  university_ed = mean(university_ed), income_top20 = mean(income_top20), Age = mean(Age)))
pred.sens.gov <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = 0, city_over200k = mean(city_over200k), 
                                  university_ed = mean(university_ed), income_top20 = mean(income_top20), Age = mean(Age)))
pred.sens.nongov <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

# Plot by supporter
pred.mf <- matrix(nrow = 4,ncol = 5)
pred.mf[c(1,3),4] <- "Gov. supporter"
pred.mf[c(2,4),4] <- "Others"
pred.mf[1:2,5] <- "List"
pred.mf[3:4,5] <- "Direct"
for (i in 1:3) {
  pred.mf[1,i] <- as.numeric(avg.pred.gov$fit[[i]])
}
for (i in 1:3) {
  pred.mf[2,i] <- as.numeric(avg.pred.nongov$fit[[i]])
}
pred.mf[3,1] <- pred.sens.gov[[1]]
pred.mf[3,2] <- pred.sens.gov[[1]]-1.96*pred.sens.gov[[2]]
pred.mf[3,3] <- pred.sens.gov[[1]]+1.96*pred.sens.gov[[2]]
pred.mf[4,1] <- pred.sens.nongov[[1]]
pred.mf[4,2] <- pred.sens.nongov[[1]]-1.96*pred.sens.nongov[[2]]
pred.mf[4,3] <- pred.sens.nongov[[1]]+1.96*pred.sens.nongov[[2]]

pred.mf <- data.frame(pred.mf)
pred.mf$X1 <- as.numeric(levels(pred.mf$X1))[pred.mf$X1]
pred.mf$X2 <- as.numeric(levels(pred.mf$X2))[pred.mf$X2]
pred.mf$X3 <- as.numeric(levels(pred.mf$X3))[pred.mf$X3]

gov.list.gr <- ggplot(pred.mf,aes(x=X4,y=X1,colour=X5)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.25) +
  geom_point(position = position_dodge(width = 0.2),size=1.3) +
  scale_y_continuous(limits=c(0, 0.83)) +
  labs(y="",x="",colour="") +
  theme_minimal(base_size = 12) + scale_colour_grey(start=0.5, end=0.2) +
  theme(axis.text.x = element_text(size=6.5), axis.title.y = element_text(size=5.5, vjust = 1), 
        axis.text.y = element_text(size=5.3), legend.text = element_text(size=6.5),
        panel.grid.minor = element_line(size = 0.15), panel.grid.major = element_line(size = 0.3))

# List analysis, Income
data_complete_notop20 <- data_complete_top20 <- data_complete
data_complete_top20[, "income_top20"] <- 1
data_complete_notop20[, "income_top20"] <- 0

avg.pred.top20 <- predict(fit.list,
                        newdata = data_complete_top20, avg = TRUE, interval = "confidence")

avg.pred.notop20 <- predict(fit.list,
                           newdata = data_complete_notop20, avg = TRUE, interval = "confidence")

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = mean(gov_support_01), city_over200k = mean(city_over200k), 
                                  university_ed = mean(university_ed), income_top20 = 1, Age = mean(Age)))
pred.sens.top20 <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

newdata2 <- with(data, data.frame(female = mean(female,na.rm = TRUE), gov_support_01 = mean(gov_support_01), city_over200k = mean(city_over200k), 
                                  university_ed = mean(university_ed), income_top20 = 0, Age = mean(Age)))
pred.sens.notop20 <- predict(fit.sens, newdata = newdata2, type = "response", se.fit = TRUE, avg = TRUE)

# Plot by income
pred.mf <- matrix(nrow = 4,ncol = 5)
pred.mf[c(1,3),4] <- "Top 20%"
pred.mf[c(2,4),4] <- "Not top 20%"
pred.mf[1:2,5] <- "List"
pred.mf[3:4,5] <- "Direct"
for (i in 1:3) {
  pred.mf[1,i] <- as.numeric(avg.pred.top20$fit[[i]])
}
for (i in 1:3) {
  pred.mf[2,i] <- as.numeric(avg.pred.notop20$fit[[i]])
}
pred.mf[3,1] <- pred.sens.top20[[1]]
pred.mf[3,2] <- pred.sens.top20[[1]]-1.96*pred.sens.top20[[2]]
pred.mf[3,3] <- pred.sens.top20[[1]]+1.96*pred.sens.top20[[2]]
pred.mf[4,1] <- pred.sens.notop20[[1]]
pred.mf[4,2] <- pred.sens.notop20[[1]]-1.96*pred.sens.notop20[[2]]
pred.mf[4,3] <- pred.sens.notop20[[1]]+1.96*pred.sens.notop20[[2]]

pred.mf <- data.frame(pred.mf)
pred.mf$X1 <- as.numeric(levels(pred.mf$X1))[pred.mf$X1]
pred.mf$X2 <- as.numeric(levels(pred.mf$X2))[pred.mf$X2]
pred.mf$X3 <- as.numeric(levels(pred.mf$X3))[pred.mf$X3]
pred.mf$X4 = factor(pred.mf$X4,levels(pred.mf$X4)[c(2,1)])

top20.list.gr <- ggplot(pred.mf,aes(x=X4,y=X1,colour=X5)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.25) +
  geom_point(position = position_dodge(width = 0.2),size=1.3) +
  scale_y_continuous(limits=c(0, 0.83)) +
  labs(y="",x="",colour="") +
  theme_minimal(base_size = 12) + scale_colour_grey(start=0.5, end=0.2) +
  theme(axis.text.x = element_text(size=6.5), axis.title.y = element_text(size=5.5, vjust = 1), 
        axis.text.y = element_text(size=5.3), legend.text = element_text(size=6.5),
        panel.grid.minor = element_line(size = 0.15), panel.grid.major = element_line(size = 0.3))

# Combine list graphs
list.combined <- ggarrange(gov.list.gr, gender.list.gr, top20.list.gr,
                          ncol = 3, nrow = 1, common.legend = TRUE, legend = "bottom")
list.combined <- annotate_figure(list.combined, left = text_grob("Estimated proportion",size = 7,rot = 90,hjust = 0.05,vjust = 3))
ggsave("listsubgroups.pdf",width=18,height=8, units = "cm")

#####################################################################################################################
# REGRESSION ESTIMATES
#####################################################################################################################
# TABLE 1 + FIGURE 2
# Corruption society
fit.soc <- lm_robust(corr_soc ~ gov_support_01 + corr_soc_treatment + gov_support_01:corr_soc_treatment, data=data)
summary(fit.soc)

data_nongov <- data_gov <- data_gov_treat <- data[1,]
data_gov_treat[, c("gov_support_01","corr_soc_treatment")] <- 1
data_gov[, "gov_support_01"] <- 1
data_gov[, "corr_soc_treatment"] <- 0
data_nongov[, c("gov_support_01","corr_soc_treatment")] <- 0

fit.soc.pred.gt <- predict(fit.soc, newdata = data_gov_treat, se.fit = TRUE, interval = "confidence", avg = TRUE)
fit.soc.pred.g <- predict(fit.soc, newdata = data_gov, se.fit = TRUE, interval = "confidence", avg = TRUE)
fit.soc.pred.ng <- predict(fit.soc, newdata = data_nongov, se.fit = TRUE, interval = "confidence", avg = TRUE)

pred.mf <- matrix(nrow = 3,ncol = 3)
pred.mf[1,] <- fit.soc.pred.ng[[1]]
pred.mf[2,] <- fit.soc.pred.g[[1]]
pred.mf[3,] <- fit.soc.pred.gt[[1]]
pred.mf <- data.frame(pred.mf)
pred.mf$est <- c("Others \n (control)","Gov. supporters \n (control)","Gov. supporters \n (prime)")
pred.mf$est <- factor(pred.mf$est)
pred.mf$est <-  factor(pred.mf$est,levels(pred.mf$est)[c(3,1,2)])

soc.gr <- ggplot(pred.mf,aes(x=est,y=X1,colour=est)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.4) +
  geom_point(position = position_dodge(width = 0.2),size=2) +
  scale_y_continuous(limits=c(1,5)) +
  labs(title="Corruption Increase During Past Year",x="",colour="",y="") +
  theme_minimal(base_size = 9) + scale_colour_grey(start=0.5, end=0.2) + 
  theme(legend.position="none", axis.text.x = element_text(size=7), axis.title.y = element_text(size=6, vjust = 1), 
        axis.text.y = element_text(size=5.5), plot.title = element_text(size=10))

# Economy
fit.econ <- lm_robust(econ ~ gov_support_01 + econ_treatment + gov_support_01:econ_treatment, data=data)
summary(fit.econ)

data_nongov <- data_gov <- data_gov_treat <- data[1,]
data_gov_treat[, c("gov_support_01","econ_treatment")] <- 1
data_gov[, "gov_support_01"] <- 1
data_gov[, "econ_treatment"] <- 0
data_nongov[, c("gov_support_01","econ_treatment")] <- 0

fit.econ.pred.gt <- predict(fit.econ, newdata = data_gov_treat, se.fit = TRUE, interval = "confidence", avg = TRUE)
fit.econ.pred.g <- predict(fit.econ, newdata = data_gov, se.fit = TRUE, interval = "confidence", avg = TRUE)
fit.econ.pred.ng <- predict(fit.econ, newdata = data_nongov, se.fit = TRUE, interval = "confidence", avg = TRUE)

pred.mf <- matrix(nrow = 3,ncol = 3)
pred.mf[1,] <- fit.econ.pred.ng[[1]]
pred.mf[2,] <- fit.econ.pred.g[[1]]
pred.mf[3,] <- fit.econ.pred.gt[[1]]
pred.mf <- data.frame(pred.mf)
pred.mf$est <- c("Others \n (control)","Gov. supporters \n (control)","Gov. supporters \n (prime)")
pred.mf$est <- factor(pred.mf$est)
pred.mf$est <-  factor(pred.mf$est,levels(pred.mf$est)[c(3,1,2)])

econ.gr <- ggplot(pred.mf,aes(x=est,y=X1,colour=est)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.4) +
  geom_point(position = position_dodge(width = 0.2),size=2) +
  scale_y_continuous(limits=c(1,5)) +
  labs(title="Worse Economy During Past Year",x="",colour="",y="") +
  theme_minimal(base_size = 9) + scale_colour_grey(start=0.5, end=0.2) + theme(legend.position="none") +
  theme(legend.position="none", axis.text.x = element_text(size=7), axis.title.y = element_text(size=6, vjust = 1), 
        axis.text.y = element_text(size=5.5), plot.title = element_text(size=10))

# Corruption in politics
fit.pol <- lm_robust(corr_pol ~ gov_support_01 + corr_pol_treatment + gov_support_01:corr_pol_treatment, data=data)
summary(fit.pol)

data_nongov <- data_gov <- data_gov_treat <- data[1,]
data_gov_treat[, c("gov_support_01","corr_pol_treatment")] <- 1
data_gov[, "gov_support_01"] <- 1
data_gov[, "corr_pol_treatment"] <- 0
data_nongov[, c("gov_support_01","corr_pol_treatment")] <- 0

fit.pol.pred.gt <- predict(fit.pol, newdata = data_gov_treat, se.fit = TRUE, interval = "confidence", avg = TRUE)
fit.pol.pred.g <- predict(fit.pol, newdata = data_gov, se.fit = TRUE, interval = "confidence", avg = TRUE)
fit.pol.pred.ng <- predict(fit.pol, newdata = data_nongov, se.fit = TRUE, interval = "confidence", avg = TRUE)

pred.mf <- matrix(nrow = 3,ncol = 3)
pred.mf[1,] <- fit.pol.pred.ng[[1]]
pred.mf[2,] <- fit.pol.pred.g[[1]]
pred.mf[3,] <- fit.pol.pred.gt[[1]]
pred.mf <- data.frame(pred.mf)
pred.mf$est <- c("Others \n (control)","Gov. supporters \n (control)","Gov. supporters \n (prime)")
pred.mf$est <- factor(pred.mf$est)
pred.mf$est <-  factor(pred.mf$est,levels(pred.mf$est)[c(3,1,2)])

pol.gr <- ggplot(pred.mf,aes(x=est,y=X1,colour=est)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.4) +
  geom_point(position = position_dodge(width = 0.2),size=2) +
  scale_y_continuous(limits=c(1,5)) +
  labs(title="Corruption in Politics",x="",colour="",y="") +
  theme_minimal(base_size = 9) + scale_colour_grey(start=0.5, end=0.2) + theme(legend.position="none") +
  theme(legend.position="none", axis.text.x = element_text(size=7), axis.title.y = element_text(size=6, vjust = 1), 
        axis.text.y = element_text(size=5.5), plot.title = element_text(size=10))

# Corruption worry
fit.worry <- lm_robust(corr_worry ~ gov_support_01 + corr_worry_treatment + gov_support_01:corr_worry_treatment, data=data)
summary(fit.worry)

data_nongov <- data_gov <- data_gov_treat <- data[1,]
data_gov_treat[, c("gov_support_01","corr_worry_treatment")] <- 1
data_gov[, "gov_support_01"] <- 1
data_gov[, "corr_worry_treatment"] <- 0
data_nongov[, c("gov_support_01","corr_worry_treatment")] <- 0

fit.worry.pred.gt <- predict(fit.worry, newdata = data_gov_treat, se.fit = TRUE, interval = "confidence", avg = TRUE)
fit.worry.pred.g <- predict(fit.worry, newdata = data_gov, se.fit = TRUE, interval = "confidence", avg = TRUE)
fit.worry.pred.ng <- predict(fit.worry, newdata = data_nongov, se.fit = TRUE, interval = "confidence", avg = TRUE)

pred.mf <- matrix(nrow = 3,ncol = 3)
pred.mf[1,] <- fit.worry.pred.ng[[1]]
pred.mf[2,] <- fit.worry.pred.g[[1]]
pred.mf[3,] <- fit.worry.pred.gt[[1]]
pred.mf <- data.frame(pred.mf)
pred.mf$est <- c("Others \n (control)","Gov. supporters \n (control)","Gov. supporters \n (prime)")
pred.mf$est <- factor(pred.mf$est)
pred.mf$est <-  factor(pred.mf$est,levels(pred.mf$est)[c(3,1,2)])

worry.gr <- ggplot(pred.mf,aes(x=est,y=X1,colour=est)) +
  geom_errorbar(aes(ymin=X2, ymax=X3), width=.1,position = position_dodge(width = 0.2), size = 0.4) +
  geom_point(position = position_dodge(width = 0.2),size=2) +
  scale_y_continuous(limits=c(1,4)) +
  labs(title="Worry About Corruption",x="",colour="", y="") +
  theme_minimal(base_size = 9) + scale_colour_grey(start=0.5, end=0.2) + theme(legend.position="none") +
  theme(legend.position="none", axis.text.x = element_text(size=7), axis.title.y = element_text(size=6, vjust = 1), 
        axis.text.y = element_text(size=5.5), plot.title = element_text(size=10))

# Combine regression graphs
reg.combined <- ggarrange(soc.gr, econ.gr, pol.gr, worry.gr,
                    ncol = 2, nrow = 2)
reg.combined <- annotate_figure(reg.combined, left = text_grob("Estimate",size = 8,rot = 90))
ggsave("regestimates.pdf")

# Regression table
fit.soc2 <- lm(corr_soc ~ gov_support_01 + corr_soc_treatment + gov_support_01:corr_soc_treatment, data=data)
fit.econ2 <- lm(econ ~ gov_support_01 + econ_treatment + gov_support_01:econ_treatment, data=data)
fit.pol2 <- lm(corr_pol ~ gov_support_01 + corr_pol_treatment + gov_support_01:corr_pol_treatment, data=data)
fit.worry2 <- lm(corr_worry ~ gov_support_01 + corr_worry_treatment + gov_support_01:corr_worry_treatment, data=data)

names(fit.soc2$coefficients)[names(fit.soc2$coefficients) == "gov_support_01"] <- "Government support"
names(fit.econ2$coefficients)[names(fit.econ2$coefficients) == "gov_support_01"] <- "Government support"
names(fit.pol2$coefficients)[names(fit.pol2$coefficients) == "gov_support_01"] <- "Government support"
names(fit.worry2$coefficients)[names(fit.worry2$coefficients) == "gov_support_01"] <- "Government support"

names(fit.soc2$coefficients)[names(fit.soc2$coefficients) == "corr_soc_treatment"] <- "Prime"
names(fit.econ2$coefficients)[names(fit.econ2$coefficients) == "econ_treatment"] <- "Prime"
names(fit.pol2$coefficients)[names(fit.pol2$coefficients) == "corr_pol_treatment"] <- "Prime"
names(fit.worry2$coefficients)[names(fit.worry2$coefficients) == "corr_worry_treatment"] <- "Prime"

names(fit.soc2$coefficients)[names(fit.soc2$coefficients) == "gov_support_01:corr_soc_treatment"] <- "Gov. support x Prime"
names(fit.econ2$coefficients)[names(fit.econ2$coefficients) == "gov_support_01:econ_treatment"] <- "Gov. support x Prime"
names(fit.pol2$coefficients)[names(fit.pol2$coefficients) == "gov_support_01:corr_pol_treatment"] <- "Gov. support x Prime"
names(fit.worry2$coefficients)[names(fit.worry2$coefficients) == "gov_support_01:corr_worry_treatment"] <- "Gov. support x Prime"

# Tablenotes were added using "threeparttable" in LaTeX
stargazer(fit.soc2,fit.econ2,fit.pol2,fit.worry2, title="The effect of the political prime on perceived corruption",
          align=TRUE, dep.var.labels.include = TRUE, 
          dep.var.labels = c("Corruption increase","Worse economy","Corruption in politics","Corruption worry"),
          omit.stat=c("LL","ser","f","rsq"), no.space=TRUE, se = starprep(fit.soc2,fit.econ2,fit.pol2,fit.worry2),
          notes = "All models are estimated using OLS. Robust standard errors in parentheses (HC2). `Prime' is an indicator variable equal to 1 if a respondent was asked the political questions before a specific corruption/economy question, and 0 otherwise. $^{*}$p$<$0.05; $^{**}$p$<$0.01; $^{***}$p$<$0.001.",  
          notes.append = FALSE, notes.align = "l", column.sep.width = "-2pt", digits = 2,
          star.cutoffs = c(0.05,0.01,0.001), font.size = "footnotesize")

